Diffusion#
import numpy as np
import scipy
import ipywidgets as widgets
import matplotlib.pyplot as plt
%matplotlib inline
Microscopic aspects of Diffusion#
Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement of particle in d dimensions to time \((t = n \Delta t)\). Generate N random walks with \(n\) steps all starting from origin \(r_0=0\).
Expressing steps in terms of time increments \(n=\frac{t}{\delta t}\) we compute average over N random walkers in \(d=1,2,3\) dimenion.
Grouping constants together we define the diffusion coefficient which is expressed in terms of microscopic quantities defined in the random walk model!
We end up with a general expression for a mean square displacement as a function of time. Any motion which adheres to this scaling with time will be called diffusive.
d=3
t= np.linspace(0, 1, 10000)
for D in [0.01, 0.1, 1, 10]:
plt.plot(t, 2*d *D * t**0.5)
#plt.loglog(t, 2*d *D * t**0.5)
plt.ylabel('$MSD(t)$')
plt.xlabel('$t$')
Text(0.5, 0, '$t$')
Brownian motion#
The Brownian motion describes the motion of a particle suspended in a solvent consisting of much smaller moleules. The displacement of particle is generatd by a sum of large number of independent random collisions with solvent molecules. Hence we invoke Central Limit Theorem to approximate displacement at each \(dt\) as a normalally idstributed random variable:
We assume we have started at position \(\mu=0\), and our variance is given by \(\sigma^2=2Dt\), Where D is the diffusion coefficient, which is related to the parameters of the discrete random walk as shown in the lecture.
In the last step, we re-wrote Brownian motion equation in a convenient way by shifting the normally distributed random variable by mean and scaling by standard deviation \(N(\mu, \sigma^2) = \mu + \sigma N(0,1)\)
def brown(T, N, dt=1, D=1, d=3):
"""
Creates 3D brownian path given:
time T
N=1 trajecotires
dt=1 timestep
D=1 diffusion coeff
returns np.array with shape (N, T, 3)
"""
n = int(T/dt) # how many points to sample
dR = np.sqrt(2*D*dt) * np.random.randn(N, n, d) # 3D position of brownian particle
R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
return R
R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
def brownian_plot(t=10):
fig, ax = plt.subplots(ncols=2)
ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
ax[1].hist(R[:, 10, 0], density=True, color='red')
ax[1].hist(R[:, t, 0], density=True)
ax[1].set_ylim([0,0.1])
ax[0].set_ylim([-200, 200])
ax[0].set_xlim([-200, 200])
fig.tight_layout()
import holoviews as hv
hv.extension('plotly')
plots = []
for i in range(10):
plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
plots.append(plot)
hv.Overlay(plots)
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]
xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])
hv.Overlay(rw_curve) << ydist << xdist
Diffusion Equation#
The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\). Formulated empirically as Fick’s laws
This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour
This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.
Can solve numerically by writing derivatives as finite differences!
Can also simulate via random walk!
Diffusion coefficient \(D\), Units \([L^2]/[T]\)
Important special case solution (here written in 1d):
where \(\sigma(t)=\sqrt{2{D}t}\)
density remains Gaussian and Gaussian becomes wider with time1
check that this is indeed a solution by plugging into the diffusion equation!
def sigma(t, D = 1):
return np.sqrt(2*D*t)
def gaussian(x, t):
return 1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
def diffusion(t=0.001):
R = brown(T=101, N=1000)
x = np.linspace(-20, 20, 100)
plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
plt.legend()
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
plt.xlim([-25, 25])